function [FixedStatesMap,AllStatesMap,SingleStateIdx] = StatesMapping(ExoFixedStates,ExoTimeStates,AgeState,SRDStates,opt)
% The indices follow a nested structure. 
% Output SingleStateIdx will give the index for all state variables considered for each observation.
% How to read that value? Columns of AllStateMap(index,:) will inform the indices, respectively, for Age, Exogenous
% Time State variables and Exogenous Fixed State Variables.
%
% Index for Exogenous Fixed State Variable (ind_f): The columns of FixedStatesMap(ind_f,:)
% informs the indices for the all exogenous fixed state variables. Those indices can be then
% mapped back into values for the state variables using
% MidPointsFixedStates, which is generated by GenExoFixedStates.m
%
% Index for Exogenous Time State Variable (ind_f): This can be mapped back
% into actual price data using pproc.y{r}, where r is the region.


% First create the indices for the exogenous fixed states:
Ns = size(ExoFixedStates,2)-1;
Ny = opt.NumFixedStates;

FixedStatesMap = MultIdx(Ny,Ns);
FixedStatesMap = [repmat(FixedStatesMap,[3 1]) kron((1:3)',ones(size(FixedStatesMap,1),1))]; 

ind = ones(size(ExoFixedStates,1),1);
for s = 1:Ns
    ind = ind + (ExoFixedStates(:,end-s)-1).*( Ny^( Ns - s ) );
end
ind = ind + (ExoFixedStates(:,end) - 1)*(Ny^Ns);

FixedExoStateIdx = ind; % Fixed Exogenous State vector index

Na = opt.MaxAge;
Np = opt.pproc.Ns^2;
Nf = size(FixedStatesMap,1);
Nd = length(opt.SRDthresholds);

MV = NaN(Na*Np*Nf*Nd,3);
tempa = (1:Na)';
tempp = (1:Np)';
tempf = (1:Nf)';
tempd = (1:Nd)';
MV(:,1) = repmat(tempa,[Np*Nf*Nd,1]); % First column with age indices
MV(:,2) = repmat(kron(tempp,ones(Na,1)),[Nf*Nd,1]); % Second column exogenous time states indices
MV(:,3) = repmat(kron(tempf,ones(Na*Np,1)),[Nd,1]); % Third column exogenous fixed states indices
MV(:,4) = repmat(kron(tempd,ones(Na*Np*Nf,1)),[1,1]); % Forth column sugarcane distance states indices

% Now classify observations in the data:
Nt = size(AgeState,2);

SingleStateIdx = NaN(size(AgeState));

for t = 1:Nt
    SingleStateIdx(:,t) = 1 + (SRDStates(:,t) - 1).*(Np*Na*Nf) + (FixedExoStateIdx - 1).*(Np*Na) + (ExoTimeStates(:,t) - 1).*(Na) + (AgeState(:,t) - 1);
end

AllStatesMap = MV;



    function MIdx = MultIdx(Ns,Ny)
        Nstar = Ns^Ny;
        MIdx = NaN(Nstar,Ny);
        UVndx = (1:Ns)';
        MIdx(:,1) = repmat(UVndx, Ns^(Ny-1), 1);

        for j = 2 : Ny
           n1          = Ns ^ (j - 1);
           n2          = Ns ^ (Ny - j);
           MIdx(:,j)  = repmat(kron(UVndx, ones(n1,1)), n2, 1);
        end
    end

end